Static eye-body allometry in Sergestidae
Katie Thomas
9/21/2021
Eye-body allometry in each species
We first examine static post-metamorphosis eye diameter-body length allometry across 16 species of deep-sea sergestids. Allometric relationships were fit using standardised major axis regression via the smatr v.3.4.8 package. For each species, standard model checks, model fits, and an interactive plot of data and fits are shown. Hover over the plots to see individual data values. Points are colored by the source of the sample (fresh collection or preserved in ethanol or paraformaldehyde).
Note that for 3 of the species, a large/unreasonable outlier was present. We identified these by eye and they have been removed from the data subsets and analyses, but the outlier value is noted under each species header for reference, and a second plot including the outlier as a red diamond is included.
#import specimen data
specimens <- data.frame(read.csv("../Data/cleaned data/specimen_data_tidy.csv", header=TRUE))#set colors for fresh and preserved samples
col_pres <- c("ethanol" = "#a6cee3",
"fresh" = "#b2df8a",
"paraformaldehyde" = "#1f78b4")Allosergestes pectinatus
# Subset data -----
apec <- specimens %>% filter(genus_species == "Allosergestes_pectinatus")
# Fit SMA model ------
sma_apec <- sma(formula = Eye_Diameter ~ Body_Length,
data = apec,
log = "xy", #sets both x and y variables as logged
method="SMA", #defines SMA as model
alpha = 0.05)
#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_apec, which = "default",type = "o")
plot(sma_apec, which = "residual",type = "o")
plot(sma_apec, which = "qq", type = "o")
par(mfrow = c(1,1)) 
#view fits
summary(sma_apec)## Call: sma(formula = Eye_Diameter ~ Body_Length, data = apec, log = "xy",
## method = "SMA", alpha = 0.05)
##
## Fit using Standardized Major Axis
##
## These variables were log-transformed before fitting: xy
##
## Confidence intervals (CI) are at 95%
##
## ------------------------------------------------------------
## Coefficients:
## elevation slope
## estimate -2.812398 1.940053
## lower limit -3.491769 1.485663
## upper limit -2.133027 2.533418
##
## H0 : variables uncorrelated
## R-squared : 0.8156794
## P-value : 9.643e-06#save coefficients of fits as object
cc_sma_apec <- data.frame(coef(sma_apec))
# Make plot ------
plot_apec <- ggplot(apec, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) +
geom_point(size = 2, alpha = 0.8) +
scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
ggtitle("Allosergestes pectinatus") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
geom_abline(data = cc_sma_apec, aes(intercept = cc_sma_apec[1,1], slope = cc_sma_apec[2,1])) +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_apec[2,1], digits = 2), sep = " "), size = 3, col = "black") +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_apec[1,1], digits = 2), sep = " "), size = 3, col = "black")
ggplotly(plot_apec)Allosergestes sargassi
Note: this species has a big outlier that was removed for analyses and plots (eye 0.7, body 13.0).
# Subset data -----
#note: outlier removed
asar <- specimens %>%
filter(genus_species == "Allosergestes_sargassi") %>%
filter(!(Eye_Diameter == round(0.67, 1) & Body_Length == round(12.97,1)))
# Fit SMA model ------
sma_asar <- sma(formula = Eye_Diameter ~ Body_Length,
data = asar,
log = "xy", #sets both x and y variables as logged
method="SMA", #defines SMA as model
alpha = 0.05)
#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_asar, which = "default",type = "o")
plot(sma_asar, which = "residual",type = "o")
plot(sma_asar, which = "qq", type = "o")
par(mfrow = c(1,1)) 
#view fits
summary(sma_asar)## Call: sma(formula = Eye_Diameter ~ Body_Length, data = asar, log = "xy",
## method = "SMA", alpha = 0.05)
##
## Fit using Standardized Major Axis
##
## These variables were log-transformed before fitting: xy
##
## Confidence intervals (CI) are at 95%
##
## ------------------------------------------------------------
## Coefficients:
## elevation slope
## estimate -2.728530 1.844597
## lower limit -3.402995 1.419494
## upper limit -2.054065 2.397008
##
## H0 : variables uncorrelated
## R-squared : 0.4037707
## P-value : 2.3852e-05#save coefficients of fits as object
cc_sma_asar <- data.frame(coef(sma_asar))
# Make plot ------
plot_asar <- ggplot(asar, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) +
geom_point(size = 2, alpha = 0.8) +
scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
ggtitle("Allosergestes sargassi") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
geom_abline(data = cc_sma_asar, aes(intercept = cc_sma_asar[1,1], slope = cc_sma_asar[2,1])) +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_asar[2,1], digits = 2), sep = " "), size = 3, col = "black") +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_asar[1,1], digits = 2), sep = " "), size = 3, col = "black")
#plot with outlier shown
plot_asar_out <- plot_asar +
geom_point(aes(y=round(0.67, 1), x=round(12.97,1)), colour="red", pch = 18) +
ggtitle("Allosergestes sargassi + outlier")
#plot
ggplotly(plot_asar)#plot with outlier
ggplotly(plot_asar_out)Challengerosergia hansjacobi
# Subset data -----
chan <- specimens %>% filter(genus_species == "Challengerosergia_hansjacobi")
# Fit SMA model ------
sma_chan <- sma(formula = Eye_Diameter ~ Body_Length,
data = chan,
log = "xy", #sets both x and y variables as logged
method="SMA", #defines SMA as model
alpha = 0.05)
#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_chan, which = "default",type = "o")
plot(sma_chan, which = "residual",type = "o")
plot(sma_chan, which = "qq", type = "o")
par(mfrow = c(1,1)) 
#view fits
summary(sma_chan)## Call: sma(formula = Eye_Diameter ~ Body_Length, data = chan, log = "xy",
## method = "SMA", alpha = 0.05)
##
## Fit using Standardized Major Axis
##
## These variables were log-transformed before fitting: xy
##
## Confidence intervals (CI) are at 95%
##
## ------------------------------------------------------------
## Coefficients:
## elevation slope
## estimate -1.528299 1.047036
## lower limit -1.913161 0.824145
## upper limit -1.143436 1.330207
##
## H0 : variables uncorrelated
## R-squared : 0.4129122
## P-value : 3.3885e-06#save coefficients of fits as object
cc_sma_chan <- data.frame(coef(sma_chan))
# Make plot ------
plot_chan <- ggplot(chan, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) +
geom_point(size = 2, alpha = 0.8) +
scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
ggtitle("Challengerosergia hansjacobi") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
geom_abline(data = cc_sma_chan, aes(intercept = cc_sma_chan[1,1], slope = cc_sma_chan[2,1])) +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_chan[2,1], digits = 2), sep = " "), size = 3, col = "black") +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_chan[1,1], digits = 2), sep = " "), size = 3, col = "black")
ggplotly(plot_chan)Challengerosergia talismani
# Subset data -----
ctal <- specimens %>% filter(genus_species == "Challengerosergia_talismani")
# Fit SMA model ------
sma_ctal <- sma(formula = Eye_Diameter ~ Body_Length,
data = ctal,
log = "xy", #sets both x and y variables as logged
method="SMA", #defines SMA as model
alpha = 0.05)
#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_ctal, which = "default",type = "o")
plot(sma_ctal, which = "residual",type = "o")
plot(sma_ctal, which = "qq", type = "o")
par(mfrow = c(1,1)) 
#view fits
summary(sma_ctal)## Call: sma(formula = Eye_Diameter ~ Body_Length, data = ctal, log = "xy",
## method = "SMA", alpha = 0.05)
##
## Fit using Standardized Major Axis
##
## These variables were log-transformed before fitting: xy
##
## Confidence intervals (CI) are at 95%
##
## ------------------------------------------------------------
## Coefficients:
## elevation slope
## estimate -1.578063 1.0389917
## lower limit -1.840235 0.8845412
## upper limit -1.315891 1.2204109
##
## H0 : variables uncorrelated
## R-squared : 0.8120894
## P-value : 2.0471e-12#save coefficients of fits as object
cc_sma_ctal <- data.frame(coef(sma_ctal))
# Make plot ------
plot_ctal <- ggplot(ctal, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) +
geom_point(size = 2, alpha = 0.8) +
scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
ggtitle("Challengerosergia talismani") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
geom_abline(data = cc_sma_ctal, aes(intercept = cc_sma_ctal[1,1], slope = cc_sma_ctal[2,1])) +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_ctal[2,1], digits = 2), sep = " "), size = 3, col = "black") +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_ctal[1,1], digits = 2), sep = " "), size = 3, col = "black")
ggplotly(plot_ctal)Deosergestes corniculum
# Subset data -----
dcor <- specimens %>% filter(genus_species == "Deosergestes_corniculum")
# Fit SMA model ------
sma_dcor <- sma(formula = Eye_Diameter ~ Body_Length,
data = dcor,
log = "xy", #sets both x and y variables as logged
method="SMA", #defines SMA as model
alpha = 0.05)
#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_dcor, which = "default",type = "o")
plot(sma_dcor, which = "residual",type = "o")
plot(sma_dcor, which = "qq", type = "o")
par(mfrow = c(1,1)) 
#view fits
summary(sma_dcor)## Call: sma(formula = Eye_Diameter ~ Body_Length, data = dcor, log = "xy",
## method = "SMA", alpha = 0.05)
##
## Fit using Standardized Major Axis
##
## These variables were log-transformed before fitting: xy
##
## Confidence intervals (CI) are at 95%
##
## ------------------------------------------------------------
## Coefficients:
## elevation slope
## estimate -1.432046 0.8653057
## lower limit -1.737583 0.7020666
## upper limit -1.126509 1.0664999
##
## H0 : variables uncorrelated
## R-squared : 0.8421211
## P-value : 8.1768e-08#save coefficients of fits as object
cc_sma_dcor <- data.frame(coef(sma_dcor))
# Make plot ------
plot_dcor <- ggplot(dcor, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) +
geom_point(size = 2, alpha = 0.8) +
scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
ggtitle("Deosergestes corniculum") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
geom_abline(data = cc_sma_dcor, aes(intercept = cc_sma_dcor[1,1], slope = cc_sma_dcor[2,1])) +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_dcor[2,1], digits = 2), sep = " "), size = 3, col = "black") +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_dcor[1,1], digits = 2), sep = " "), size = 3, col = "black")
ggplotly(plot_dcor)Deosergestes henseni
# Subset data -----
dhen <- specimens %>% filter(genus_species == "Deosergestes_henseni")
# Fit SMA model ------
sma_dhen <- sma(formula = Eye_Diameter ~ Body_Length,
data = dhen,
log = "xy", #sets both x and y variables as logged
method="SMA", #defines SMA as model
alpha = 0.05)
#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_dhen, which = "default",type = "o")
plot(sma_dhen, which = "residual",type = "o")
plot(sma_dhen, which = "qq", type = "o")
par(mfrow = c(1,1)) 
#view fits
summary(sma_dhen)## Call: sma(formula = Eye_Diameter ~ Body_Length, data = dhen, log = "xy",
## method = "SMA", alpha = 0.05)
##
## Fit using Standardized Major Axis
##
## These variables were log-transformed before fitting: xy
##
## Confidence intervals (CI) are at 95%
##
## ------------------------------------------------------------
## Coefficients:
## elevation slope
## estimate -1.306994 0.8238705
## lower limit -1.593941 0.6637531
## upper limit -1.020047 1.0226131
##
## H0 : variables uncorrelated
## R-squared : 0.3380975
## P-value : 1.6982e-06#save coefficients of fits as object
cc_sma_dhen <- data.frame(coef(sma_dhen))
# Make plot ------
plot_dhen <- ggplot(dhen, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) +
geom_point(size = 2, alpha = 0.8) +
scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
ggtitle("Deosergestes henseni") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
geom_abline(data = cc_sma_dhen, aes(intercept = cc_sma_dhen[1,1], slope = cc_sma_dhen[2,1])) +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_dhen[2,1], digits = 2), sep = " "), size = 3, col = "black") +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_dhen[1,1], digits = 2), sep = " "), size = 3, col = "black")
ggplotly(plot_dhen)Eusergestes arcticus
# Subset data -----
euar <- specimens %>% filter(genus_species == "Eusergestes_arcticus")
# Fit SMA model ------
sma_euar <- sma(formula = Eye_Diameter ~ Body_Length,
data = euar,
log = "xy", #sets both x and y variables as logged
method="SMA", #defines SMA as model
alpha = 0.05)
#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_euar, which = "default",type = "o")
plot(sma_euar, which = "residual",type = "o")
plot(sma_euar, which = "qq", type = "o")
par(mfrow = c(1,1)) 
#view fits
summary(sma_euar)## Call: sma(formula = Eye_Diameter ~ Body_Length, data = euar, log = "xy",
## method = "SMA", alpha = 0.05)
##
## Fit using Standardized Major Axis
##
## These variables were log-transformed before fitting: xy
##
## Confidence intervals (CI) are at 95%
##
## ------------------------------------------------------------
## Coefficients:
## elevation slope
## estimate -2.1539322 1.3306114
## lower limit -3.8635100 0.6694485
## upper limit -0.4443543 2.6447543
##
## H0 : variables uncorrelated
## R-squared : 0.8368056
## P-value : 0.029484#save coefficients of fits as object
cc_sma_euar <- data.frame(coef(sma_euar))
# Make plot ------
plot_euar <- ggplot(euar, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) +
geom_point(size = 2, alpha = 0.8) +
scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
ggtitle("Eusergestes arcticus") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
geom_abline(data = cc_sma_euar, aes(intercept = cc_sma_euar[1,1], slope = cc_sma_euar[2,1])) +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_euar[2,1], digits = 2), sep = " "), size = 3, col = "black") +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_euar[1,1], digits = 2), sep = " "), size = 3, col = "black")
ggplotly(plot_euar)Gardinerosergia splendens
Note: this species has a big outlier that was removed for analyses and plots (eye 1.7, body 22.0).
# Subset data -----
gspl <- specimens %>%
filter(genus_species == "Gardinerosergia_splendens") %>%
filter(!(Eye_Diameter == round(1.71, 1) & Body_Length == round(21.98,1)))
# Fit SMA model ------
sma_gspl <- sma(formula = Eye_Diameter ~ Body_Length,
data = gspl,
log = "xy", #sets both x and y variables as logged
method="SMA", #defines SMA as model
alpha = 0.05)
#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_gspl, which = "default",type = "o")
plot(sma_gspl, which = "residual",type = "o")
plot(sma_gspl, which = "qq", type = "o")
par(mfrow = c(1,1)) 
#view fits
summary(sma_gspl)## Call: sma(formula = Eye_Diameter ~ Body_Length, data = gspl, log = "xy",
## method = "SMA", alpha = 0.05)
##
## Fit using Standardized Major Axis
##
## These variables were log-transformed before fitting: xy
##
## Confidence intervals (CI) are at 95%
##
## ------------------------------------------------------------
## Coefficients:
## elevation slope
## estimate -1.2253950 0.8523703
## lower limit -1.5326713 0.6732695
## upper limit -0.9181186 1.0791149
##
## H0 : variables uncorrelated
## R-squared : 0.4597308
## P-value : 1.1132e-06#save coefficients of fits as object
cc_sma_gspl <- data.frame(coef(sma_gspl))
# Make plot ------
plot_gspl <- ggplot(gspl, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) +
geom_point(size = 2, alpha = 0.8) +
scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
ggtitle("Gardinerosergia splendens") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
geom_abline(data = cc_sma_gspl, aes(intercept = cc_sma_gspl[1,1], slope = cc_sma_gspl[2,1])) +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_gspl[2,1], digits = 2), sep = " "), size = 3, col = "black") +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_gspl[1,1], digits = 2), sep = " "), size = 3, col = "black")
#plot with outlier shown
plot_gspl_out <- plot_gspl +
geom_point(aes(y=round(1.71,1), x=round(21.98,1)), colour="red", pch = 18) +
ggtitle("Gardinerosergia splendens + outlier")
#plot
ggplotly(plot_gspl)#plot with outlier
ggplotly(plot_gspl_out)Neosergestes edwardsii
Note: this species has a big outlier that was removed for analyses and plots (eye 0.9, body 0.7). Also note that for this species, sample size was insufficient to estimate eye-body allometry (N = 6, p = 0.57 for correlation between eye size and body size).
# Subset data -----
nedw <- specimens %>%
filter(genus_species == "Neosergestes_edwardsii") %>%
filter(!(Eye_Diameter == round(0.85,1) & Body_Length == round(0.72,1)))
# Fit SMA model ------
sma_nedw <- sma(formula = Eye_Diameter ~ Body_Length,
data = nedw,
log = "xy", #sets both x and y variables as logged
method="SMA", #defines SMA as model
alpha = 0.05)
#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_nedw, which = "default",type = "o")
plot(sma_nedw, which = "residual",type = "o")
plot(sma_nedw, which = "qq", type = "o")
par(mfrow = c(1,1)) 
#view fits
summary(sma_nedw)## Call: sma(formula = Eye_Diameter ~ Body_Length, data = nedw, log = "xy",
## method = "SMA", alpha = 0.05)
##
## Fit using Standardized Major Axis
##
## These variables were log-transformed before fitting: xy
##
## Confidence intervals (CI) are at 95%
##
## ------------------------------------------------------------
## Coefficients:
## elevation slope
## estimate 1.4244363 -1.2712830
## lower limit -0.7533312 -3.8007834
## upper limit 3.6022037 -0.4252177
##
## H0 : variables uncorrelated
## R-squared : 0.08540156
## P-value : 0.57413#save coefficients of fits as object
cc_sma_nedw <- data.frame(coef(sma_nedw))
# Make plot ------
plot_nedw <- ggplot(nedw, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) +
geom_point(size = 2, alpha = 0.8) +
scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
ggtitle("Neosergestes edwardsii") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic"))
#plot with outlier shown
plot_nedw_out <- plot_nedw +
geom_point(aes(x=round(0.85,1), y=round(0.72,1)), colour="red", pch = 18) +
ggtitle("Neosergestes edwardsii + outlier")
#plot
ggplotly(plot_nedw)#plot with outlier
ggplotly(plot_nedw_out)Parasergestes armatus
# Subset data -----
parm <- specimens %>% filter(genus_species == "Parasergestes_armatus")
# Fit SMA model ------
sma_parm <- sma(formula = Eye_Diameter ~ Body_Length,
data = parm,
log = "xy", #sets both x and y variables as logged
method="SMA", #defines SMA as model
alpha = 0.05)
#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_parm, which = "default",type = "o")
plot(sma_parm, which = "residual",type = "o")
plot(sma_parm, which = "qq", type = "o")
par(mfrow = c(1,1)) 
#view fits
summary(sma_parm)## Call: sma(formula = Eye_Diameter ~ Body_Length, data = parm, log = "xy",
## method = "SMA", alpha = 0.05)
##
## Fit using Standardized Major Axis
##
## These variables were log-transformed before fitting: xy
##
## Confidence intervals (CI) are at 95%
##
## ------------------------------------------------------------
## Coefficients:
## elevation slope
## estimate -1.697577 1.0573117
## lower limit -2.095490 0.8217616
## upper limit -1.299663 1.3603801
##
## H0 : variables uncorrelated
## R-squared : 0.5670657
## P-value : 1.5744e-06#save coefficients of fits as object
cc_sma_parm <- data.frame(coef(sma_parm))
# Make plot ------
plot_parm <- ggplot(parm, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) +
geom_point(size = 2, alpha = 0.8) +
scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
ggtitle("Parasergestes armatus") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
geom_abline(data = cc_sma_parm, aes(intercept = cc_sma_parm[1,1], slope = cc_sma_parm[2,1])) +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_parm[2,1], digits = 2), sep = " "), size = 3, col = "black") +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_parm[1,1], digits = 2), sep = " "), size = 3, col = "black")
ggplotly(plot_parm)Parasergestes vigilax
# Subset data -----
pvig <- specimens %>% filter(genus_species == "Parasergestes_vigilax")
# Fit SMA model ------
sma_pvig <- sma(formula = Eye_Diameter ~ Body_Length,
data = pvig,
log = "xy", #sets both x and y variables as logged
method="SMA", #defines SMA as model
alpha = 0.05)
#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_pvig, which = "default",type = "o")
plot(sma_pvig, which = "residual",type = "o")
plot(sma_pvig, which = "qq", type = "o")
par(mfrow = c(1,1)) 
#view fits
summary(sma_pvig)## Call: sma(formula = Eye_Diameter ~ Body_Length, data = pvig, log = "xy",
## method = "SMA", alpha = 0.05)
##
## Fit using Standardized Major Axis
##
## These variables were log-transformed before fitting: xy
##
## Confidence intervals (CI) are at 95%
##
## ------------------------------------------------------------
## Coefficients:
## elevation slope
## estimate -3.306489 2.375713
## lower limit -4.455793 1.641298
## upper limit -2.157186 3.438750
##
## H0 : variables uncorrelated
## R-squared : 0.5274941
## P-value : 0.00096145#save coefficients of fits as object
cc_sma_pvig <- data.frame(coef(sma_pvig))
# Make plot ------
plot_pvig <- ggplot(pvig, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) +
geom_point(size = 2, alpha = 0.8) +
scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
ggtitle("Parasergestes vigilax") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
geom_abline(data = cc_sma_pvig, aes(intercept = cc_sma_pvig[1,1], slope = cc_sma_pvig[2,1])) +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_pvig[2,1], digits = 2), sep = " "), size = 3, col = "black") +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_pvig[1,1], digits = 2), sep = " "), size = 3, col = "black")
ggplotly(plot_pvig)Phorcosergia grandis
# Subset data -----
pgra <- specimens %>% filter(genus_species == "Phorcosergia_grandis")
# Fit SMA model ------
sma_pgra <- sma(formula = Eye_Diameter ~ Body_Length,
data = pgra,
log = "xy", #sets both x and y variables as logged
method="SMA", #defines SMA as model
alpha = 0.05)
#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_pgra, which = "default",type = "o")
plot(sma_pgra, which = "residual",type = "o")
plot(sma_pgra, which = "qq", type = "o")
par(mfrow = c(1,1)) 
#view fits
summary(sma_pgra)## Call: sma(formula = Eye_Diameter ~ Body_Length, data = pgra, log = "xy",
## method = "SMA", alpha = 0.05)
##
## Fit using Standardized Major Axis
##
## These variables were log-transformed before fitting: xy
##
## Confidence intervals (CI) are at 95%
##
## ------------------------------------------------------------
## Coefficients:
## elevation slope
## estimate -1.470665 0.9366565
## lower limit -1.632022 0.8480626
## upper limit -1.309308 1.0345055
##
## H0 : variables uncorrelated
## R-squared : 0.8798127
## P-value : < 2.22e-16#save coefficients of fits as object
cc_sma_pgra <- data.frame(coef(sma_pgra))
# Make plot ------
plot_pgra <- ggplot(pgra, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) +
geom_point(size = 2, alpha = 0.8) +
scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
ggtitle("Phorcosergia grandis") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
geom_abline(data = cc_sma_pgra, aes(intercept = cc_sma_pgra[1,1], slope = cc_sma_pgra[2,1])) +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_pgra[2,1], digits = 2), sep = " "), size = 3, col = "black") +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_pgra[1,1], digits = 2), sep = " "), size = 3, col = "black")
ggplotly(plot_pgra)Robustosergia robusta
# Subset data -----
rrob <- specimens %>% filter(genus_species == "Robustosergia_robusta")
# Fit SMA model ------
sma_rrob <- sma(formula = Eye_Diameter ~ Body_Length,
data = rrob,
log = "xy", #sets both x and y variables as logged
method="SMA", #defines SMA as model
alpha = 0.05)
#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_rrob, which = "default",type = "o")
plot(sma_rrob, which = "residual",type = "o")
plot(sma_rrob, which = "qq", type = "o")
par(mfrow = c(1,1)) 
#view fits
summary(sma_rrob)## Call: sma(formula = Eye_Diameter ~ Body_Length, data = rrob, log = "xy",
## method = "SMA", alpha = 0.05)
##
## Fit using Standardized Major Axis
##
## These variables were log-transformed before fitting: xy
##
## Confidence intervals (CI) are at 95%
##
## ------------------------------------------------------------
## Coefficients:
## elevation slope
## estimate -1.2073093 0.8386026
## lower limit -1.5732188 0.6553974
## upper limit -0.8413998 1.0730198
##
## H0 : variables uncorrelated
## R-squared : 0.6989428
## P-value : 6.7833e-07#save coefficients of fits as object
cc_sma_rrob <- data.frame(coef(sma_rrob))
# Make plot ------
plot_rrob <- ggplot(rrob, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) +
geom_point(size = 2, alpha = 0.8) +
scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
ggtitle("Robustosergia robusta") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
geom_abline(data = cc_sma_rrob, aes(intercept = cc_sma_rrob[1,1], slope = cc_sma_rrob[2,1])) +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_rrob[2,1], digits = 2), sep = " "), size = 3, col = "black") +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_rrob[1,1], digits = 2), sep = " "), size = 3, col = "black")
ggplotly(plot_rrob)Robustosergia regalis
# Subset data -----
rreg <- specimens %>% filter(genus_species == "Robustosergia_regalis")
# Fit SMA model ------
sma_rreg <- sma(formula = Eye_Diameter ~ Body_Length,
data = rreg,
log = "xy", #sets both x and y variables as logged
method="SMA", #defines SMA as model
alpha = 0.05)
#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_rreg, which = "default",type = "o")
plot(sma_rreg, which = "residual",type = "o")
plot(sma_rreg, which = "qq", type = "o")
par(mfrow = c(1,1)) 
#view fits
summary(sma_rreg)## Call: sma(formula = Eye_Diameter ~ Body_Length, data = rreg, log = "xy",
## method = "SMA", alpha = 0.05)
##
## Fit using Standardized Major Axis
##
## These variables were log-transformed before fitting: xy
##
## Confidence intervals (CI) are at 95%
##
## ------------------------------------------------------------
## Coefficients:
## elevation slope
## estimate -1.324440 0.8978007
## lower limit -1.531510 0.7807317
## upper limit -1.117369 1.0324239
##
## H0 : variables uncorrelated
## R-squared : 0.8280327
## P-value : 2.4971e-15#save coefficients of fits as object
cc_sma_rreg <- data.frame(coef(sma_rreg))
# Make plot ------
plot_rreg <- ggplot(rreg, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) +
geom_point(size = 2, alpha = 0.8) +
scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
ggtitle("Robustosergia regalis") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
geom_abline(data = cc_sma_rreg, aes(intercept = cc_sma_rreg[1,1], slope = cc_sma_rreg[2,1])) +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_rreg[2,1], digits = 2), sep = " "), size = 3, col = "black") +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_rreg[1,1], digits = 2), sep = " "), size = 3, col = "black")
ggplotly(plot_rreg)Sergestes atlanticus
# Subset data -----
satl <- specimens %>% filter(genus_species == "Sergestes_atlanticus")
# Fit SMA model ------
sma_satl <- sma(formula = Eye_Diameter ~ Body_Length,
data = satl,
log = "xy", #sets both x and y variables as logged
method="SMA", #defines SMA as model
alpha = 0.05)
#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_satl, which = "default",type = "o")
plot(sma_satl, which = "residual",type = "o")
plot(sma_satl, which = "qq", type = "o")
par(mfrow = c(1,1)) 
#view fits
summary(sma_satl)## Call: sma(formula = Eye_Diameter ~ Body_Length, data = satl, log = "xy",
## method = "SMA", alpha = 0.05)
##
## Fit using Standardized Major Axis
##
## These variables were log-transformed before fitting: xy
##
## Confidence intervals (CI) are at 95%
##
## ------------------------------------------------------------
## Coefficients:
## elevation slope
## estimate -2.839127 1.951572
## lower limit -3.727822 1.418514
## upper limit -1.950431 2.684945
##
## H0 : variables uncorrelated
## R-squared : 0.7067609
## P-value : 8.6609e-05#save coefficients of fits as object
cc_sma_satl <- data.frame(coef(sma_satl))
# Make plot ------
plot_satl <- ggplot(satl, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) +
geom_point(size = 2, alpha = 0.8) +
scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
ggtitle("Sergestes atlanticus") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
geom_abline(data = cc_sma_satl, aes(intercept = cc_sma_satl[1,1], slope = cc_sma_satl[2,1])) +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_satl[2,1], digits = 2), sep = " "), size = 3, col = "black") +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_satl[1,1], digits = 2), sep = " "), size = 3, col = "black")
ggplotly(plot_satl)Sergia tenuiremis
# Subset data -----
sten <- specimens %>% filter(genus_species == "Sergia_tenuiremis")
# Fit SMA model ------
sma_sten <- sma(formula = Eye_Diameter ~ Body_Length,
data = sten,
log = "xy", #sets both x and y variables as logged
method="SMA", #defines SMA as model
alpha = 0.05)
#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_sten, which = "default",type = "o")
plot(sma_sten, which = "residual",type = "o")
plot(sma_sten, which = "qq", type = "o")
par(mfrow = c(1,1)) 
#view fits
summary(sma_sten)## Call: sma(formula = Eye_Diameter ~ Body_Length, data = sten, log = "xy",
## method = "SMA", alpha = 0.05)
##
## Fit using Standardized Major Axis
##
## These variables were log-transformed before fitting: xy
##
## Confidence intervals (CI) are at 95%
##
## ------------------------------------------------------------
## Coefficients:
## elevation slope
## estimate -1.2313110 0.8025570
## lower limit -1.6822103 0.5747311
## upper limit -0.7804117 1.1206941
##
## H0 : variables uncorrelated
## R-squared : 0.2880646
## P-value : 0.0032339#save coefficients of fits as object
cc_sma_sten <- data.frame(coef(sma_sten))
# Make plot ------
plot_sten <- ggplot(sten, aes(x = Body_Length, y = Eye_Diameter, col = Preservation)) +
geom_point(size = 2, alpha = 0.8) +
scale_color_manual(values = col_pres, name = "Sample state", breaks = c("ethanol", "fresh", "paraformaldehyde")) +
scale_x_log10(name = "Body length (mm)") + #makes x axis log scale and named
scale_y_log10(name = "Eye diameter (mm)") + #makes y axis log scale and named
ggtitle("Sergia tenuiremis") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(face = "italic")) +
geom_abline(data = cc_sma_sten, aes(intercept = cc_sma_sten[1,1], slope = cc_sma_sten[2,1])) +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 1.5, nudge_x = 4, label = paste("slope =", round(cc_sma_sten[2,1], digits = 2), sep = " "), size = 3, col = "black") +
geom_text(x = -Inf, y = Inf, hjust = -0.05, vjust = 3, label = paste("intercept =", round(cc_sma_sten[1,1], digits = 2), sep = " "), size = 3, col = "black")
ggplotly(plot_sten)Figure for static eye-body scaling across sergestid species
#make figure panels
fig.a <- plot_dcor + theme(legend.position = "none", axis.title.x = element_blank())
fig.b <- plot_dhen + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
fig.c <- plot_apec + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
fig.d <- plot_asar + theme(legend.position = "none", axis.title.x = element_blank())
fig.e <- plot_satl + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
fig.f <- plot_pvig + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
fig.g <- plot_parm + theme(legend.position = "none", axis.title.x = element_blank())
fig.h <- plot_euar + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
fig.i <- plot_gspl + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
fig.j <- plot_rreg + theme(legend.position = "none", axis.title.x = element_blank())
fig.k <- plot_rrob + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
fig.l <- plot_pgra + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
fig.m <- plot_sten + theme(legend.position = "none")
fig.n <- plot_ctal + theme(legend.position = "none", axis.title.y = element_blank())
fig.o <- plot_chan + theme(legend.position = "none", axis.title.y = element_blank())
# arrange panels in figure
plots <- plot_grid(fig.a, fig.b, fig.c, fig.d, fig.e, fig.f, fig.g, fig.h, fig.i, fig.j, fig.k, fig.l, fig.m, fig.n, fig.o, #list of plots to arrange in grid
align = 'vh', #align horizontally and vertically
axis = 'lb',
labels = c("A", "B", "C", "D", "E", "F", "G", "H", "I","J", "K", "L", "M", "N", "O"), #panel labels for figure
hjust = -1, #adjustment for panel labels
nrow = 5) #number of rows in grids
#extract legend
leg <- get_legend(fig.e + theme(legend.position="bottom"))
#view figure
plot_grid(plots, leg, nrow = 2, rel_heights = c(1, .01), align = 'vh')
#export figure
#png("../Figures/dev-allometry.png", width = 10, height = 12, units = "in", res = 500)
pdf("../Figures/Fig-S1.pdf", width = 10, height = 12)
plot_grid(plots, leg, nrow = 2, rel_heights = c(1, .05), align = 'vh')
dev.off()Test whether species differ in static allometry
We’ve looked at how eyes scale with body size within each species, but are there significant differences in static allometry across species? Here we use two approaches: first, we compare linear mixed models to see which best describes our data. Then, we look at pairwise comparisons of SMA regressions among species to see which are driving these differences.
Linear mixed model to test whether slopes vary across species
Here, we use linear mixed models in lme4 v.1.1.25 to test whether species show differences in static eye-body allometric slopes. We fit 1) a constant-slope, variable intercept model to our eye and body size data with species as a random effect and 2) a variable slope model with the same effects, and then compared model fits by AIC scores. Both models are fit using reduced maximum likelihood (REML), as this approach is better than maximum likelihood (ML) for detecting differences in random effects among models, which is what we are testing here (ML is better for fixed effects).
First we removed the three outliers and the species with insufficient sampling (N. edwardsii) from our full specimen dataset, then ran each model.
#Remove outliers from specimen dataset
specimens_ed <- specimens %>%
filter(!(genus_species == "Allosergestes_sargassi" & Eye_Diameter == 0.7 & Body_Length == 13.0)) %>% #(0.67, 12.97 raw)
filter(!(genus_species == "Gardinerosergia_splendens" & Eye_Diameter == 1.7 & Body_Length == 22.0)) %>% #(1.71, 21.98 raw)
filter(!(genus_species == "Neosergestes_edwardsii" & Eye_Diameter == 0.8 & Body_Length == 0.7)) #(0.85, 0.72 raw)
#Export data without outliers
write.csv(specimens_ed, file = "../Data/cleaned data/specimen_data_tidy_no-outliers.csv", row.names = FALSE)
#Remove species with insufficient sampling
specimens_test <- specimens_ed %>%
filter(genus_species != "Neosergestes_edwardsii") %>%
mutate(species_text = as.factor(paste(Genus, Species, sep = " "))) # Add labels for species textFixed slope (variable intercepts) model
# variable intercepts model (species can have different intercepts but share mean slope) -------
lme_fixed <- lmer(log10(Eye_Diameter) ~ log10(Body_Length) + (1|genus_species),
data = specimens_test)
# summary
summary(lme_fixed)## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(Eye_Diameter) ~ log10(Body_Length) + (1 | genus_species)
## Data: specimens_test
##
## REML criterion at convergence: -1091.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.6547 -0.5660 0.0425 0.5699 3.7587
##
## Random effects:
## Groups Name Variance Std.Dev.
## genus_species (Intercept) 0.004257 0.06525
## Residual 0.004549 0.06745
## Number of obs: 450, groups: genus_species, 15
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -1.24683 0.05121 -24.35
## log10(Body_Length) 0.80642 0.03122 25.83
##
## Correlation of Fixed Effects:
## (Intr)
## lg10(Bdy_L) -0.941# intercepts (variable) and slope (constant)
coef(lme_fixed)## $genus_species
## (Intercept) log10(Body_Length)
## Allosergestes_pectinatus -1.337988 0.8064169
## Allosergestes_sargassi -1.295962 0.8064169
## Challengerosergia_hansjacobi -1.164606 0.8064169
## Challengerosergia_talismani -1.216875 0.8064169
## Deosergestes_corniculum -1.328655 0.8064169
## Deosergestes_henseni -1.278582 0.8064169
## Eusergestes_arcticus -1.246866 0.8064169
## Gardinerosergia_splendens -1.158266 0.8064169
## Parasergestes_armatus -1.325352 0.8064169
## Parasergestes_vigilax -1.299386 0.8064169
## Phorcosergia_grandis -1.246064 0.8064169
## Robustosergia_regalis -1.176361 0.8064169
## Robustosergia_robusta -1.155267 0.8064169
## Sergestes_atlanticus -1.234190 0.8064169
## Sergia_tenuiremis -1.238014 0.8064169
##
## attr(,"class")
## [1] "coef.mer"# AIC
AIC(lme_fixed)## [1] -1083.888#fixed effects
fixef(lme_fixed)## (Intercept) log10(Body_Length)
## -1.2468288 0.8064169#confidence intervals
## if these overlap zero, we do not have good support for significant effects
## if .sig01, which is our estimate of the variability in the random effect is very large and very widely defined, t indicates we may have a lack of precision between our groups - either because the group effect is small between groups, we have too few groups to get a more precise estimate, we have too few units within each group, or a combination of all of the above.
confint(lme_fixed, level = 0.95)## 2.5 % 97.5 %
## .sig01 0.04428477 0.09578174
## .sigma 0.06313279 0.07211335
## (Intercept) -1.34880013 -1.14643360
## log10(Body_Length) 0.74536494 0.86929119Variable slopes model
# variable slopes model (species can have different slopes/interaction effect between body size and species on eye size) -------
lme_variable <- lmer(log10(Eye_Diameter) ~ log10(Body_Length) + (1 + log10(Body_Length) | genus_species),
data = specimens_test)
# summary
summary(lme_variable)## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(Eye_Diameter) ~ log10(Body_Length) + (1 + log10(Body_Length) |
## genus_species)
## Data: specimens_test
##
## REML criterion at convergence: -1125.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.4279 -0.5359 0.0306 0.5910 4.2687
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## genus_species (Intercept) 0.232807 0.48250
## log10(Body_Length) 0.104532 0.32331 -0.99
## Residual 0.003977 0.06306
## Number of obs: 450, groups: genus_species, 15
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -1.35693 0.13739 -9.876
## log10(Body_Length) 0.89778 0.09181 9.778
##
## Correlation of Fixed Effects:
## (Intr)
## lg10(Bdy_L) -0.993# intercepts (variable) and slope (constant)
coef(lme_variable)## $genus_species
## (Intercept) log10(Body_Length)
## Allosergestes_pectinatus -2.2120964 1.4762026
## Allosergestes_sargassi -1.7513868 1.1366814
## Challengerosergia_hansjacobi -1.0760213 0.7484245
## Challengerosergia_talismani -1.4111004 0.9318641
## Deosergestes_corniculum -1.1953553 0.7273404
## Deosergestes_henseni -0.7843247 0.4963146
## Eusergestes_arcticus -1.2644263 0.8196651
## Gardinerosergia_splendens -0.8716955 0.6169987
## Parasergestes_armatus -1.3258605 0.8075505
## Parasergestes_vigilax -2.0269469 1.3719045
## Phorcosergia_grandis -1.3536721 0.8693323
## Robustosergia_regalis -1.2328227 0.8412590
## Robustosergia_robusta -1.1498050 0.8037740
## Sergestes_atlanticus -1.8945693 1.2766884
## Sergia_tenuiremis -0.8038697 0.5427274
##
## attr(,"class")
## [1] "coef.mer"# AIC
AIC(lme_variable)## [1] -1113.879#fixed effects
fixef(lme_variable)## (Intercept) log10(Body_Length)
## -1.3569302 0.8977818#confidence intervals
## if these overlap zero, we do not have good support for significant effects
## if .sig01, which is our estimate of the variability in the random effect is very large and very widely defined, t indicates we may have a lack of precision between our groups - either because the group effect is small between groups, we have too few groups to get a more precise estimate, we have too few units within each group, or a combination of all of the above.
confint(lme_variable, level = 0.95)## 2.5 % 97.5 %
## .sig01 0.27582256 0.74974111
## .sig02 -0.99580186 -0.97348803
## .sig03 0.17701121 0.50987554
## .sigma 0.05903897 0.06767603
## (Intercept) -1.64344107 -1.07983292
## log10(Body_Length) 0.71455352 1.09409513Note that the parameter estimates for slopes and intercepts for species groups are a bit different than those we estimated via OLS. That’s because here, the model is looking across all the data for all the species, and finding a common mean slope/intercept as well as the group effects.
Compare models
Then we can compare the two models via AIC and log likelihood.
#model AIC comparison
AIC(lme_fixed, lme_variable)## df AIC
## lme_fixed 4 -1083.888
## lme_variable 6 -1113.879#log likelihood comparison
anova(lme_fixed, lme_variable, refit=FALSE)## Data: specimens_test
## Models:
## lme_fixed: log10(Eye_Diameter) ~ log10(Body_Length) + (1 | genus_species)
## lme_variable: log10(Eye_Diameter) ~ log10(Body_Length) + (1 + log10(Body_Length) | genus_species)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## lme_fixed 4 -1083.9 -1067.5 545.94 -1091.9
## lme_variable 6 -1113.9 -1089.2 562.94 -1125.9 33.991 2 4.159e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Here we find by comparing AIC scores and log likelihoods that the model with variable slopes is better supported than the one with fixed slopes. We have evidence here that species differ significantly in the slopes of their static eye-body allometries rather than following a shared growth trajectory.
Test for significant effect of preservation as a covariate
Here, we add preservation method as a fixed effect into our linear mixed models of eye size vs. body size with species as a random effect. We have already selected the variable slopes model, so here we use a maximum likelihood approach (ML) to fit models that 1) don’t include preservation as a covariate and 2) do include preservation as a covariate and then compare model fits. We include preservation as a fixed rather than a random effect because random effects should have more than 5 levels, and preservation type has only 3.
Model without preservation type
# variable slopes model (species can have different slopes/interaction effect between body size and species on eye size) -------
lme_nopres <- lmer(log10(Eye_Diameter) ~ log10(Body_Length) + (1 + log10(Body_Length) | genus_species),
data = specimens_test,
REML = FALSE) #uses ML instead of REML
# summary
summary(lme_nopres)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: log10(Eye_Diameter) ~ log10(Body_Length) + (1 + log10(Body_Length) |
## genus_species)
## Data: specimens_test
##
## AIC BIC logLik deviance df.resid
## -1123.3 -1098.6 567.6 -1135.3 444
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.4474 -0.5397 0.0290 0.5871 4.2627
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## genus_species (Intercept) 0.205885 0.45375
## log10(Body_Length) 0.091871 0.30310 -0.99
## Residual 0.003983 0.06311
## Number of obs: 450, groups: genus_species, 15
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -1.35193 0.13036 -10.37
## log10(Body_Length) 0.89371 0.08686 10.29
##
## Correlation of Fixed Effects:
## (Intr)
## lg10(Bdy_L) -0.993# intercepts (variable) and slope (constant)
coef(lme_nopres)## $genus_species
## (Intercept) log10(Body_Length)
## Allosergestes_pectinatus -2.1804353 1.4518078
## Allosergestes_sargassi -1.7432761 1.1308509
## Challengerosergia_hansjacobi -1.0817350 0.7520768
## Challengerosergia_talismani -1.4101465 0.9312293
## Deosergestes_corniculum -1.1915651 0.7253703
## Deosergestes_henseni -0.7877696 0.4985142
## Eusergestes_arcticus -1.2643317 0.8197744
## Gardinerosergia_splendens -0.8765164 0.6200697
## Parasergestes_armatus -1.3286935 0.8097102
## Parasergestes_vigilax -1.9939834 1.3460569
## Phorcosergia_grandis -1.3527442 0.8688288
## Robustosergia_regalis -1.2356712 0.8429113
## Robustosergia_robusta -1.1602404 0.8095839
## Sergestes_atlanticus -1.8546790 1.2481074
## Sergia_tenuiremis -0.8172119 0.5507995
##
## attr(,"class")
## [1] "coef.mer"# AIC
AIC(lme_nopres)## [1] -1123.278#fixed effects
fixef(lme_nopres)## (Intercept) log10(Body_Length)
## -1.3519333 0.8937128#confidence intervals
## if these overlap zero, we do not have good support for significant effects
## if .sig01, which is our estimate of the variability in the random effect is very large and very widely defined, t indicates we may have a lack of precision between our groups - either because the group effect is small between groups, we have too few groups to get a more precise estimate, we have too few units within each group, or a combination of all of the above.
confint(lme_nopres, level = 0.95)## 2.5 % 97.5 %
## .sig01 0.27582268 0.74974134
## .sig02 -0.99580165 -0.97348804
## .sig03 0.17701128 0.50987568
## .sigma 0.05903897 0.06767603
## (Intercept) -1.64344107 -1.07983293
## log10(Body_Length) 0.71455335 1.09409513Model including preservation type
#reorder levels for preservation so that contrasts will be compared to fresh samples
specimens_test$Preservation <- factor(specimens_test$Preservation,
levels = c("fresh", "ethanol", "paraformaldehyde"))
# variable slopes model (species can have different slopes/interaction effect between body size and species on eye size) -------
lme_withpres <- lmer(log10(Eye_Diameter) ~ log10(Body_Length) + Preservation + (1 + log10(Body_Length) | genus_species),
data = specimens_test,
REML = FALSE) #uses ML instead of REML
# summary
summary(lme_withpres)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: log10(Eye_Diameter) ~ log10(Body_Length) + Preservation + (1 +
## log10(Body_Length) | genus_species)
## Data: specimens_test
##
## AIC BIC logLik deviance df.resid
## -1152.9 -1120.0 584.5 -1168.9 442
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.7037 -0.5122 0.0295 0.5644 4.3676
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## genus_species (Intercept) 0.167876 0.40973
## log10(Body_Length) 0.071632 0.26764 -0.99
## Residual 0.003701 0.06083
## Number of obs: 450, groups: genus_species, 15
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -1.353719 0.119794 -11.300
## log10(Body_Length) 0.894613 0.078169 11.445
## Preservationethanol -0.019607 0.008151 -2.405
## Preservationparaformaldehyde 0.045520 0.012167 3.741
##
## Correlation of Fixed Effects:
## (Intr) l10(B_ Prsrvtnt
## lg10(Bdy_L) -0.990
## Prsrvtnthnl -0.057 0.022
## Prsrvtnprfr -0.068 0.041 0.455# AIC
AIC(lme_withpres)## [1] -1152.918#fixed effects
fixef(lme_withpres)## (Intercept) log10(Body_Length)
## -1.35371863 0.89461257
## Preservationethanol Preservationparaformaldehyde
## -0.01960702 0.04551984#confidence intervals
## if these overlap zero, we do not have good support for significant effects
## if .sig01, which is our estimate of the variability in the random effect is very large and very widely defined, t indicates we may have a lack of precision between our groups - either because the group effect is small between groups, we have too few groups to get a more precise estimate, we have too few units within each group, or a combination of all of the above.
confint(lme_withpres, level = 0.95)## 2.5 % 97.5 %
## .sig01 0.25281957 0.67323848
## .sig02 -0.99341603 -0.96585932
## .sig03 0.15905144 0.44860338
## .sigma 0.05692248 0.06521126
## (Intercept) -1.61782966 -1.10480110
## log10(Body_Length) 0.73401477 1.07261071
## Preservationethanol -0.03573800 -0.00350345
## Preservationparaformaldehyde 0.02141721 0.06962976Compare models
Then we can compare the two models via AIC and log likelihood.
#model AIC comparison
AIC(lme_nopres, lme_withpres)## df AIC
## lme_nopres 6 -1123.278
## lme_withpres 8 -1152.918#log likelihood comparison
anova(lme_nopres, lme_withpres)## Data: specimens_test
## Models:
## lme_nopres: log10(Eye_Diameter) ~ log10(Body_Length) + (1 + log10(Body_Length) | genus_species)
## lme_withpres: log10(Eye_Diameter) ~ log10(Body_Length) + Preservation + (1 + log10(Body_Length) | genus_species)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## lme_nopres 6 -1123.3 -1098.6 567.64 -1135.3
## lme_withpres 8 -1152.9 -1120.0 584.46 -1168.9 33.64 2 4.957e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Here we see that the model that includes preservation type as a covariate is significantly better supported. Ethanol preservation has a negative effect on intercept (eye size compared to body size), but paraformaldehyde has a positive effect. Note also that these effects are very small, so while they are significantly affecting our models, they seem to have minor effects.
Linear SMA regression for pairwise comparisons among species slopes
The smatr package also allows you to test whether groups exhibit significant differences in allometric scaling (slopes or intercepts). Below, we again use standardized major axis regression to model the scaling of eye diameter with body length across all specimens measured and test whether we see significant differences in static allometry across species.
# Model for static allometry across specimens with species and species*body size interaction as covariates
sma_species <- sma(formula = Eye_Diameter ~ Body_Length * genus_species,
data = specimens_test,
log = "xy", #sets both x and y variables as logged
method="SMA", #defines SMA as model
multcomp = TRUE, # adds pairwise comparisons between groups
multcompmethod = "adjusted", # adjusts p-value for multiple comparisons
alpha = 0.05)
#plot fit, residuals, qq
par(mfrow = c(2,2)) #make plot window 2x2
plot(sma_species, which = "default",type = "o")
plot(sma_species,which = "residual",type = "o")
plot(sma_species, which = "qq", type = "o")
par(mfrow = c(1,1)) 
#view fits
summary(sma_species)## Call: sma(formula = Eye_Diameter ~ Body_Length * genus_species, data = specimens_test,
## log = "xy", method = "SMA", alpha = 0.05, multcomp = TRUE,
## multcompmethod = "adjusted")
##
## Fit using Standardized Major Axis
##
## These variables were log-transformed before fitting: xy
##
## Confidence intervals (CI) are at 95%
##
## ------------------------------------------------------------
## Results of comparing lines among groups.
##
## H0 : slopes are equal.
## Likelihood ratio statistic : 78.37 with 14 degrees of freedom
## P-value : 5.6738e-11
## ------------------------------------------------------------
##
## Results of multiple comparisons among groups.
##
## Test for pair-wise difference in slope :
## genus_species_1 genus_species_2 Pval
## 1 Allosergestes_pectinatus Allosergestes_sargassi 1.00000000
## 2 Allosergestes_pectinatus Challengerosergia_hansjacobi 0.10631965
## 3 Allosergestes_pectinatus Challengerosergia_talismani 0.03771725
## 4 Allosergestes_pectinatus Deosergestes_corniculum 0.00370375
## 5 Allosergestes_pectinatus Deosergestes_henseni 0.00191135
## 6 Allosergestes_pectinatus Eusergestes_arcticus 1.00000000
## 7 Allosergestes_pectinatus Gardinerosergia_splendens 0.00417512
## 8 Allosergestes_pectinatus Parasergestes_armatus 0.14662731
## 9 Allosergestes_pectinatus Parasergestes_vigilax 1.00000000
## 10 Allosergestes_pectinatus Phorcosergia_grandis 0.00484593
## 11 Allosergestes_pectinatus Robustosergia_regalis 0.00330606
## 12 Allosergestes_pectinatus Robustosergia_robusta 0.00386781
## 13 Allosergestes_pectinatus Sergestes_atlanticus 1.00000000
## 14 Allosergestes_pectinatus Sergia_tenuiremis 0.00977338
## 15 Allosergestes_sargassi Challengerosergia_hansjacobi 0.17642048
## 16 Allosergestes_sargassi Challengerosergia_talismani 0.03323729
## 17 Allosergestes_sargassi Deosergestes_corniculum 0.00219121
## 18 Allosergestes_sargassi Deosergestes_henseni 0.00062265
## 19 Allosergestes_sargassi Eusergestes_arcticus 1.00000000
## 20 Allosergestes_sargassi Gardinerosergia_splendens 0.00279183
## 21 Allosergestes_sargassi Parasergestes_armatus 0.25381226
## 22 Allosergestes_sargassi Parasergestes_vigilax 1.00000000
## 23 Allosergestes_sargassi Phorcosergia_grandis 0.00071873
## 24 Allosergestes_sargassi Robustosergia_regalis 0.00055799
## 25 Allosergestes_sargassi Robustosergia_robusta 0.00340588
## 26 Allosergestes_sargassi Sergestes_atlanticus 1.00000000
## 27 Allosergestes_sargassi Sergia_tenuiremis 0.01682182
## 28 Challengerosergia_hansjacobi Challengerosergia_talismani 1.00000000
## 29 Challengerosergia_hansjacobi Deosergestes_corniculum 1.00000000
## 30 Challengerosergia_hansjacobi Deosergestes_henseni 0.99999989
## 31 Challengerosergia_hansjacobi Eusergestes_arcticus 1.00000000
## 32 Challengerosergia_hansjacobi Gardinerosergia_splendens 1.00000000
## 33 Challengerosergia_hansjacobi Parasergestes_armatus 1.00000000
## 34 Challengerosergia_hansjacobi Parasergestes_vigilax 0.04921452
## 35 Challengerosergia_hansjacobi Phorcosergia_grandis 1.00000000
## 36 Challengerosergia_hansjacobi Robustosergia_regalis 1.00000000
## 37 Challengerosergia_hansjacobi Robustosergia_robusta 1.00000000
## 38 Challengerosergia_hansjacobi Sergestes_atlanticus 0.25203409
## 39 Challengerosergia_hansjacobi Sergia_tenuiremis 1.00000000
## 40 Challengerosergia_talismani Deosergestes_corniculum 0.99999999
## 41 Challengerosergia_talismani Deosergestes_henseni 0.99993757
## 42 Challengerosergia_talismani Eusergestes_arcticus 1.00000000
## 43 Challengerosergia_talismani Gardinerosergia_splendens 1.00000000
## 44 Challengerosergia_talismani Parasergestes_armatus 1.00000000
## 45 Challengerosergia_talismani Parasergestes_vigilax 0.02124209
## 46 Challengerosergia_talismani Phorcosergia_grandis 1.00000000
## 47 Challengerosergia_talismani Robustosergia_regalis 1.00000000
## 48 Challengerosergia_talismani Robustosergia_robusta 0.99999992
## 49 Challengerosergia_talismani Sergestes_atlanticus 0.11697344
## 50 Challengerosergia_talismani Sergia_tenuiremis 0.99999999
## 51 Deosergestes_corniculum Deosergestes_henseni 1.00000000
## 52 Deosergestes_corniculum Eusergestes_arcticus 0.99999999
## 53 Deosergestes_corniculum Gardinerosergia_splendens 1.00000000
## 54 Deosergestes_corniculum Parasergestes_armatus 1.00000000
## 55 Deosergestes_corniculum Parasergestes_vigilax 0.00236873
## 56 Deosergestes_corniculum Phorcosergia_grandis 1.00000000
## 57 Deosergestes_corniculum Robustosergia_regalis 1.00000000
## 58 Deosergestes_corniculum Robustosergia_robusta 1.00000000
## 59 Deosergestes_corniculum Sergestes_atlanticus 0.01350778
## 60 Deosergestes_corniculum Sergia_tenuiremis 1.00000000
## 61 Deosergestes_henseni Eusergestes_arcticus 0.99999959
## 62 Deosergestes_henseni Gardinerosergia_splendens 1.00000000
## 63 Deosergestes_henseni Parasergestes_armatus 0.99999978
## 64 Deosergestes_henseni Parasergestes_vigilax 0.00128436
## 65 Deosergestes_henseni Phorcosergia_grandis 1.00000000
## 66 Deosergestes_henseni Robustosergia_regalis 1.00000000
## 67 Deosergestes_henseni Robustosergia_robusta 1.00000000
## 68 Deosergestes_henseni Sergestes_atlanticus 0.00731827
## 69 Deosergestes_henseni Sergia_tenuiremis 1.00000000
## 70 Eusergestes_arcticus Gardinerosergia_splendens 0.99999997
## 71 Eusergestes_arcticus Parasergestes_armatus 1.00000000
## 72 Eusergestes_arcticus Parasergestes_vigilax 0.99996770
## 73 Eusergestes_arcticus Phorcosergia_grandis 1.00000000
## 74 Eusergestes_arcticus Robustosergia_regalis 1.00000000
## 75 Eusergestes_arcticus Robustosergia_robusta 0.99999993
## 76 Eusergestes_arcticus Sergestes_atlanticus 1.00000000
## 77 Eusergestes_arcticus Sergia_tenuiremis 0.99999911
## 78 Gardinerosergia_splendens Parasergestes_armatus 1.00000000
## 79 Gardinerosergia_splendens Parasergestes_vigilax 0.00266011
## 80 Gardinerosergia_splendens Phorcosergia_grandis 1.00000000
## 81 Gardinerosergia_splendens Robustosergia_regalis 1.00000000
## 82 Gardinerosergia_splendens Robustosergia_robusta 1.00000000
## 83 Gardinerosergia_splendens Sergestes_atlanticus 0.01499340
## 84 Gardinerosergia_splendens Sergia_tenuiremis 1.00000000
## 85 Parasergestes_armatus Parasergestes_vigilax 0.06386567
## 86 Parasergestes_armatus Phorcosergia_grandis 1.00000000
## 87 Parasergestes_armatus Robustosergia_regalis 1.00000000
## 88 Parasergestes_armatus Robustosergia_robusta 1.00000000
## 89 Parasergestes_armatus Sergestes_atlanticus 0.31517941
## 90 Parasergestes_armatus Sergia_tenuiremis 1.00000000
## 91 Parasergestes_vigilax Phorcosergia_grandis 0.00320670
## 92 Parasergestes_vigilax Robustosergia_regalis 0.00220257
## 93 Parasergestes_vigilax Robustosergia_robusta 0.00240837
## 94 Parasergestes_vigilax Sergestes_atlanticus 1.00000000
## 95 Parasergestes_vigilax Sergia_tenuiremis 0.00481337
## 96 Phorcosergia_grandis Robustosergia_regalis 1.00000000
## 97 Phorcosergia_grandis Robustosergia_robusta 1.00000000
## 98 Phorcosergia_grandis Sergestes_atlanticus 0.01851043
## 99 Phorcosergia_grandis Sergia_tenuiremis 1.00000000
## 100 Robustosergia_regalis Robustosergia_robusta 1.00000000
## 101 Robustosergia_regalis Sergestes_atlanticus 0.01271145
## 102 Robustosergia_regalis Sergia_tenuiremis 1.00000000
## 103 Robustosergia_robusta Sergestes_atlanticus 0.01358845
## 104 Robustosergia_robusta Sergia_tenuiremis 1.00000000
## 105 Sergestes_atlanticus Sergia_tenuiremis 0.02662028
## TestStat
## 1 7.666232e-02
## 2 1.070239e+01
## 3 1.269776e+01
## 4 1.710663e+01
## 5 1.836679e+01
## 6 1.514543e+00
## 7 1.687873e+01
## 8 1.006766e+01
## 9 8.341304e-01
## 10 1.659541e+01
## 11 1.732280e+01
## 12 1.702417e+01
## 13 8.875795e-04
## 14 1.526295e+01
## 15 9.695878e+00
## 16 1.293861e+01
## 17 1.810626e+01
## 18 2.051006e+01
## 19 1.218474e+00
## 20 1.764469e+01
## 21 8.943571e+00
## 22 1.268181e+00
## 23 2.023537e+01
## 24 2.072006e+01
## 25 1.726619e+01
## 26 7.817618e-02
## 27 1.423248e+01
## 28 2.850713e-03
## 29 1.459201e+00
## 30 2.159382e+00
## 31 7.165782e-01
## 32 1.474497e+00
## 33 3.154991e-03
## 34 1.218978e+01
## 35 7.281210e-01
## 36 1.215204e+00
## 37 1.670749e+00
## 38 8.958462e+00
## 39 1.652888e+00
## 40 1.971056e+00
## 41 2.909088e+00
## 42 7.908634e-01
## 43 1.901819e+00
## 44 1.381453e-02
## 45 1.378961e+01
## 46 1.203054e+00
## 47 1.877144e+00
## 48 2.130775e+00
## 49 1.051505e+01
## 50 1.917202e+00
## 51 1.091006e-01
## 52 1.922892e+00
## 53 9.383239e-03
## 54 1.528670e+00
## 55 1.795778e+01
## 56 4.919054e-01
## 57 9.043887e-02
## 58 3.949879e-02
## 59 1.464885e+01
## 60 1.486450e-01
## 61 2.285326e+00
## 62 4.454856e-02
## 63 2.224689e+00
## 64 1.912552e+01
## 65 1.138002e+00
## 66 4.397057e-01
## 67 1.181780e-02
## 68 1.581216e+01
## 69 1.724888e-02
## 70 2.036867e+00
## 71 6.530716e-01
## 72 2.808168e+00
## 73 1.445758e+00
## 74 1.716984e+00
## 75 2.123385e+00
## 76 1.501074e+00
## 77 2.362523e+00
## 78 1.546238e+00
## 79 1.773675e+01
## 80 5.367445e-01
## 81 1.427310e-01
## 82 9.239058e-03
## 83 1.445085e+01
## 84 8.641470e-02
## 85 1.169009e+01
## 86 7.995621e-01
## 87 1.285411e+00
## 88 1.740929e+00
## 89 8.475706e+00
## 90 1.723105e+00
## 91 1.738089e+01
## 92 1.809640e+01
## 93 1.792616e+01
## 94 6.837798e-01
## 95 1.660823e+01
## 96 2.447817e-01
## 97 7.055402e-01
## 98 1.405093e+01
## 99 7.801888e-01
## 100 2.376709e-01
## 101 1.476415e+01
## 102 3.827565e-01
## 103 1.463755e+01
## 104 4.529580e-02
## 105 1.336087e+01
##
## ------------------------------------------------------------
## Coefficients by group in variable "genus_species"
##
## Group: Allosergestes_pectinatus
## elevation slope
## estimate -2.812398 1.940053
## lower limit -3.491769 1.485663
## upper limit -2.133027 2.533418
##
## H0 : variables uncorrelated.
## R-squared : 0.8156794
## P-value : 9.643e-06
##
## Group: Allosergestes_sargassi
## elevation slope
## estimate -2.728530 1.844597
## lower limit -3.402995 1.419494
## upper limit -2.054065 2.397008
##
## H0 : variables uncorrelated.
## R-squared : 0.4037707
## P-value : 2.3852e-05
##
## Group: Challengerosergia_hansjacobi
## elevation slope
## estimate -1.528299 1.047036
## lower limit -1.913161 0.824145
## upper limit -1.143436 1.330207
##
## H0 : variables uncorrelated.
## R-squared : 0.4129122
## P-value : 3.3885e-06
##
## Group: Challengerosergia_talismani
## elevation slope
## estimate -1.578063 1.0389917
## lower limit -1.840235 0.8845412
## upper limit -1.315891 1.2204109
##
## H0 : variables uncorrelated.
## R-squared : 0.8120894
## P-value : 2.0471e-12
##
## Group: Deosergestes_corniculum
## elevation slope
## estimate -1.432046 0.8653057
## lower limit -1.737583 0.7020666
## upper limit -1.126509 1.0664999
##
## H0 : variables uncorrelated.
## R-squared : 0.8421211
## P-value : 8.1768e-08
##
## Group: Deosergestes_henseni
## elevation slope
## estimate -1.306994 0.8238705
## lower limit -1.593941 0.6637531
## upper limit -1.020047 1.0226131
##
## H0 : variables uncorrelated.
## R-squared : 0.3380975
## P-value : 1.6982e-06
##
## Group: Eusergestes_arcticus
## elevation slope
## estimate -2.1539322 1.3306114
## lower limit -3.8635100 0.6694485
## upper limit -0.4443543 2.6447543
##
## H0 : variables uncorrelated.
## R-squared : 0.8368056
## P-value : 0.029484
##
## Group: Gardinerosergia_splendens
## elevation slope
## estimate -1.2253950 0.8523703
## lower limit -1.5326713 0.6732695
## upper limit -0.9181186 1.0791149
##
## H0 : variables uncorrelated.
## R-squared : 0.4597308
## P-value : 1.1132e-06
##
## Group: Parasergestes_armatus
## elevation slope
## estimate -1.697577 1.0573117
## lower limit -2.095490 0.8217616
## upper limit -1.299663 1.3603801
##
## H0 : variables uncorrelated.
## R-squared : 0.5670657
## P-value : 1.5744e-06
##
## Group: Parasergestes_vigilax
## elevation slope
## estimate -3.306489 2.375713
## lower limit -4.455793 1.641298
## upper limit -2.157186 3.438750
##
## H0 : variables uncorrelated.
## R-squared : 0.5274941
## P-value : 0.00096145
##
## Group: Phorcosergia_grandis
## elevation slope
## estimate -1.470665 0.9366565
## lower limit -1.632022 0.8480626
## upper limit -1.309308 1.0345055
##
## H0 : variables uncorrelated.
## R-squared : 0.8798127
## P-value : < 2.22e-16
##
## Group: Robustosergia_regalis
## elevation slope
## estimate -1.324440 0.8978007
## lower limit -1.531510 0.7807317
## upper limit -1.117369 1.0324239
##
## H0 : variables uncorrelated.
## R-squared : 0.8280327
## P-value : 2.4971e-15
##
## Group: Robustosergia_robusta
## elevation slope
## estimate -1.2073093 0.8386026
## lower limit -1.5732188 0.6553974
## upper limit -0.8413998 1.0730198
##
## H0 : variables uncorrelated.
## R-squared : 0.6989428
## P-value : 6.7833e-07
##
## Group: Sergestes_atlanticus
## elevation slope
## estimate -2.839127 1.951572
## lower limit -3.727822 1.418514
## upper limit -1.950431 2.684945
##
## H0 : variables uncorrelated.
## R-squared : 0.7067609
## P-value : 8.6609e-05
##
## Group: Sergia_tenuiremis
## elevation slope
## estimate -1.2313110 0.8025570
## lower limit -1.6822103 0.5747311
## upper limit -0.7804117 1.1206941
##
## H0 : variables uncorrelated.
## R-squared : 0.2880646
## P-value : 0.0032339#save coefficients of fits as object
cc_species <- data.frame(coef(sma_species))Among 15 sergestid species (450 specimens), species significantly differ in eye-body allometric scaling (p < 0.0001). Pairwise comparisons show which species differ significantly in slope and which do not.
Here we plot the static allometry of all the species with sufficient sampling (n = 15). It’s interactive, so you can click on different species in the legend to hide or show those points.
#reorder factor levels for figure legend (phylo order)
specimens_test$species_text <- factor(specimens_test$species_text,
levels = c("Deosergestes corniculum",
"Deosergestes henseni",
"Allosergestes pectinatus",
"Allosergestes sargassi",
"Sergestes atlanticus",
"Parasergestes vigilax",
"Parasergestes armatus",
"Eusergestes arcticus",
"Gardinerosergia splendens",
"Robustosergia regalis",
"Robustosergia robusta",
"Phorcosergia grandis",
"Sergia tenuiremis",
"Challengerosergia talismani",
"Challengerosergia hansjacobi"))
#make shape pallette
shapes.sp <- c("Deosergestes corniculum" = 21,
"Deosergestes henseni" = 22,
"Allosergestes pectinatus" = 23,
"Allosergestes sargassi" = 24,
"Sergestes atlanticus" = 25,
"Parasergestes vigilax" = 21,
"Parasergestes armatus" = 22,
"Eusergestes arcticus" = 23,
"Gardinerosergia splendens" = 24,
"Robustosergia regalis" = 25,
"Robustosergia robusta" = 21,
"Phorcosergia grandis" = 22,
"Sergia tenuiremis" = 23,
"Challengerosergia talismani" = 24,
"Challengerosergia hansjacobi" = 25)
#rainbow palette
cols.sp <- c("Deosergestes corniculum" = "#16ABA8",
"Deosergestes henseni" = "#1f78b4",
"Allosergestes pectinatus" = "#3acf75",
"Allosergestes sargassi" = "#33a02c",
"Sergestes atlanticus" = "#fb9a99",
"Parasergestes vigilax" = "#e31a1c",
"Parasergestes armatus" = "#fdbf6f",
"Eusergestes arcticus" = "#ff7f00",
"Gardinerosergia splendens" = "#cab2d6",
"Robustosergia regalis" = "#6a3d9a",
"Robustosergia robusta" = "gray31",
"Phorcosergia grandis" = "burlywood4",
"Sergia tenuiremis" = "black",
"Challengerosergia talismani" = "maroon1",
"Challengerosergia hansjacobi" = "orchid1")
#sergia/sergestes green purple pallette
# cols.sp <- c(#Sergestes group
# "Deosergestes corniculum" = "#512E5F",
# "Deosergestes henseni" = "#633974",
# "Allosergestes pectinatus" = "#76448A",
# "Allosergestes sargassi" = "#884EA0",
# "Sergestes atlanticus" = "#9B59B6",
# #"Neosergestes edwardsii" = "#AF7AC5"
# "Parasergestes vigilax" = "#C39BD3",
# "Parasergestes armatus" = "#D7BDE2",
# "Eusergestes arcticus" = "#EBDEF0",
# #Sergia group
# "Gardinerosergia splendens" = "#D4EFDF",
# "Robustosergia regalis" = "#A9DFBF",
# "Robustosergia robusta" = "#52BE80",
# "Phorcosergia grandis" = "#27AE60",
# "Sergia tenuiremis" = "#1E8449",
# "Challengerosergia talismani" = "#196F3D",
# "Challengerosergia hansjacobi" = "#145A32")
#set s limits by species for plots
limits <- specimens_test %>%
group_by(genus_species) %>%
summarise(max = max(Body_Length, na.rm = TRUE),
min = min(Body_Length, na.rm = TRUE))
#make plot
plot_species <- ggplot(specimens_test, aes(x = Body_Length, y = Eye_Diameter, color = species_text, shape = species_text, fill = species_text)) +
geom_point(size = 2, alpha = 0.4) +
scale_shape_manual(values = shapes.sp, name = "Species") +
scale_color_manual(values = cols.sp, name = "Species") +
scale_fill_manual(values = cols.sp, name = "Species") +
scale_x_log10(name = "Body length (mm)", breaks = c(15,20,30,50,100)) + #makes x axis log scale and named
scale_y_log10(name = "Eye diameter (mm)", breaks = c(0.2,0.3,0.5,1, 2.6)) + #makes y axis log scale and named
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.text = element_text(face = "italic")) +
geom_segment(aes(x = limits$min[which(limits$genus_species == "Allosergestes_pectinatus")], xend = limits$max[which(limits$genus_species == "Allosergestes_pectinatus")], y = 10^(log10(limits$min[which(limits$genus_species == "Allosergestes_pectinatus")])*cc_species["Allosergestes_pectinatus", "slope"] + cc_species["Allosergestes_pectinatus", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Allosergestes_pectinatus")])*cc_species["Allosergestes_pectinatus", "slope"] + cc_species["Allosergestes_pectinatus", "elevation"])), color = cols.sp["Allosergestes pectinatus"]) + #Rana temporaria adult scaling
geom_segment(aes(x = limits$min[which(limits$genus_species == "Allosergestes_sargassi")], xend = limits$max[which(limits$genus_species == "Allosergestes_sargassi")], y = 10^(log10(limits$min[which(limits$genus_species == "Allosergestes_sargassi")])*cc_species["Allosergestes_sargassi", "slope"] + cc_species["Allosergestes_sargassi", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Allosergestes_sargassi")])*cc_species["Allosergestes_sargassi", "slope"] + cc_species["Allosergestes_sargassi", "elevation"])), color = cols.sp["Allosergestes sargassi"]) + #Allosergestes_sargassi adult scaling
geom_segment(aes(x = limits$min[which(limits$genus_species == "Challengerosergia_hansjacobi")], xend = limits$max[which(limits$genus_species == "Challengerosergia_hansjacobi")], y = 10^(log10(limits$min[which(limits$genus_species == "Challengerosergia_hansjacobi")])*cc_species["Challengerosergia_hansjacobi", "slope"] + cc_species["Challengerosergia_hansjacobi", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Challengerosergia_hansjacobi")])*cc_species["Challengerosergia_hansjacobi", "slope"] + cc_species["Challengerosergia_hansjacobi", "elevation"])), color = cols.sp["Challengerosergia hansjacobi"]) + #Challengerosergia_hansjacobi adult scaling
geom_segment(aes(x = limits$min[which(limits$genus_species == "Challengerosergia_talismani")], xend = limits$max[which(limits$genus_species == "Challengerosergia_talismani")], y = 10^(log10(limits$min[which(limits$genus_species == "Challengerosergia_talismani")])*cc_species["Challengerosergia_talismani", "slope"] + cc_species["Challengerosergia_talismani", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Challengerosergia_talismani")])*cc_species["Challengerosergia_talismani", "slope"] + cc_species["Challengerosergia_talismani", "elevation"])), color = cols.sp["Challengerosergia talismani"]) + #Challengerosergia_talismani adult scaling
geom_segment(aes(x = limits$min[which(limits$genus_species == "Deosergestes_corniculum")], xend = limits$max[which(limits$genus_species == "Deosergestes_corniculum")], y = 10^(log10(limits$min[which(limits$genus_species == "Deosergestes_corniculum")])*cc_species["Deosergestes_corniculum", "slope"] + cc_species["Deosergestes_corniculum", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Deosergestes_corniculum")])*cc_species["Deosergestes_corniculum", "slope"] + cc_species["Deosergestes_corniculum", "elevation"])), color = cols.sp["Deosergestes corniculum"]) + #Deosergestes_corniculum adult scaling
geom_segment(aes(x = limits$min[which(limits$genus_species == "Deosergestes_henseni")], xend = limits$max[which(limits$genus_species == "Deosergestes_henseni")], y = 10^(log10(limits$min[which(limits$genus_species == "Deosergestes_henseni")])*cc_species["Deosergestes_henseni", "slope"] + cc_species["Deosergestes_henseni", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Deosergestes_henseni")])*cc_species["Deosergestes_henseni", "slope"] + cc_species["Deosergestes_henseni", "elevation"])), color = cols.sp["Deosergestes henseni"]) + #Deosergestes henseni adult scaling
geom_segment(aes(x = limits$min[which(limits$genus_species == "Gardinerosergia_splendens")], xend = limits$max[which(limits$genus_species == "Gardinerosergia_splendens")], y = 10^(log10(limits$min[which(limits$genus_species == "Gardinerosergia_splendens")])*cc_species["Gardinerosergia_splendens", "slope"] + cc_species["Gardinerosergia_splendens", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Gardinerosergia_splendens")])*cc_species["Gardinerosergia_splendens", "slope"] + cc_species["Gardinerosergia_splendens", "elevation"])), color = cols.sp["Gardinerosergia splendens"]) + #Gardinerosergia splendens adult scaling
geom_segment(aes(x = limits$min[which(limits$genus_species == "Parasergestes_armatus")], xend = limits$max[which(limits$genus_species == "Parasergestes_armatus")], y = 10^(log10(limits$min[which(limits$genus_species == "Parasergestes_armatus")])*cc_species["Parasergestes_armatus", "slope"] + cc_species["Parasergestes_armatus", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Parasergestes_armatus")])*cc_species["Parasergestes_armatus", "slope"] + cc_species["Parasergestes_armatus", "elevation"])), color = cols.sp["Parasergestes armatus"]) + #Parasergestes armatus adult scaling
geom_segment(aes(x = limits$min[which(limits$genus_species == "Parasergestes_vigilax")], xend = limits$max[which(limits$genus_species == "Parasergestes_vigilax")], y = 10^(log10(limits$min[which(limits$genus_species == "Parasergestes_vigilax")])*cc_species["Parasergestes_vigilax", "slope"] + cc_species["Parasergestes_vigilax", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Parasergestes_vigilax")])*cc_species["Parasergestes_vigilax", "slope"] + cc_species["Parasergestes_vigilax", "elevation"])), color = cols.sp["Parasergestes vigilax"]) + #Parasergestes vigilax adult scaling
geom_segment(aes(x = limits$min[which(limits$genus_species == "Phorcosergia_grandis")], xend = limits$max[which(limits$genus_species == "Phorcosergia_grandis")], y = 10^(log10(limits$min[which(limits$genus_species == "Phorcosergia_grandis")])*cc_species["Phorcosergia_grandis", "slope"] + cc_species["Phorcosergia_grandis", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Phorcosergia_grandis")])*cc_species["Phorcosergia_grandis", "slope"] + cc_species["Phorcosergia_grandis", "elevation"])), color = cols.sp["Phorcosergia grandis"]) + #Phorcosergia grandis adult scaling
geom_segment(aes(x = limits$min[which(limits$genus_species == "Robustosergia_robusta")], xend = limits$max[which(limits$genus_species == "Robustosergia_robusta")], y = 10^(log10(limits$min[which(limits$genus_species == "Robustosergia_robusta")])*cc_species["Robustosergia_robusta", "slope"] + cc_species["Robustosergia_robusta", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Robustosergia_robusta")])*cc_species["Robustosergia_robusta", "slope"] + cc_species["Robustosergia_robusta", "elevation"])), color = cols.sp["Robustosergia robusta"]) + #Robustosergia robusta adult scaling
geom_segment(aes(x = limits$min[which(limits$genus_species == "Robustosergia_regalis")], xend = limits$max[which(limits$genus_species == "Robustosergia_regalis")], y = 10^(log10(limits$min[which(limits$genus_species == "Robustosergia_regalis")])*cc_species["Robustosergia_regalis", "slope"] + cc_species["Robustosergia_regalis", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Robustosergia_regalis")])*cc_species["Robustosergia_regalis", "slope"] + cc_species["Robustosergia_regalis", "elevation"])), color = cols.sp["Robustosergia regalis"]) +
geom_segment(aes(x = limits$min[which(limits$genus_species == "Sergestes_atlanticus")], xend = limits$max[which(limits$genus_species == "Sergestes_atlanticus")], y = 10^(log10(limits$min[which(limits$genus_species == "Sergestes_atlanticus")])*cc_species["Sergestes_atlanticus", "slope"] + cc_species["Sergestes_atlanticus", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Sergestes_atlanticus")])*cc_species["Sergestes_atlanticus", "slope"] + cc_species["Sergestes_atlanticus", "elevation"])), color = cols.sp["Sergestes atlanticus"]) +
geom_segment(aes(x = limits$min[which(limits$genus_species == "Sergia_tenuiremis")], xend = limits$max[which(limits$genus_species == "Sergia_tenuiremis")], y = 10^(log10(limits$min[which(limits$genus_species == "Sergia_tenuiremis")])*cc_species["Sergia_tenuiremis", "slope"] + cc_species["Sergia_tenuiremis", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Sergia_tenuiremis")])*cc_species["Sergia_tenuiremis", "slope"] + cc_species["Sergia_tenuiremis", "elevation"])), color = cols.sp["Sergia tenuiremis"])+
geom_segment(aes(x = limits$min[which(limits$genus_species == "Eusergestes_arcticus")], xend = limits$max[which(limits$genus_species == "Eusergestes_arcticus")], y = 10^(log10(limits$min[which(limits$genus_species == "Eusergestes_arcticus")])*cc_species["Eusergestes_arcticus", "slope"] + cc_species["Eusergestes_arcticus", "elevation"]), yend = 10^(log10(limits$max[which(limits$genus_species == "Eusergestes_arcticus")])*cc_species["Eusergestes_arcticus", "slope"] + cc_species["Eusergestes_arcticus", "elevation"])), color = cols.sp["Eusergestes arcticus"])+
theme(legend.text = element_text(face = "italic"))
#interactive plot
ggplotly(plot_species)#make plot with alpha = 1 for legend
plot_species_leg <- ggplot(specimens_test, aes(x = Body_Length, y = Eye_Diameter, color = species_text, shape = species_text, fill = species_text)) +
geom_point(size = 2, alpha = 1) +
scale_shape_manual(values = shapes.sp, name = "Species") +
scale_color_manual(values = cols.sp, name = "Species") +
scale_fill_manual(values = cols.sp, name = "Species") +
theme_bw() +
theme(legend.text = element_text(face = "italic"))
# remove legend from previous plot
plot <- plot_species + theme(legend.position = "none")
# extract legend from this plot with alpha = 1
leg <- get_legend(plot_species_leg)
#export figure
#png("../Figures/dev-species.png", width = 12, height = 12, units = "in", res = 500)
pdf("../Figures/Fig-1.pdf", width = 9, height = 5)
plot_grid(plot, leg, ncol = 2, rel_widths = c(1, .4))
dev.off()